#ifndef PBVP_TETRA_H
#define PBVP_TETRA_H

#include "../math/Vector3.h"
#include "../data/geometrydata.h"
#include "edgecache.h"

namespace PBVP {
	/// A class representing Tetrahedron
	/**
	* contains functions for handling function integral on tetrahedrons
	*/
	class Tetra
	{
	public:
		/// Constructor
		/**
		*  The vertices with be sorted by their values such that v1 < v2
		*  < v3 < v4
		*
		* @param x1 coordinates of vertex 1
		* @param x2 coordinates of vertex 2
		* @param x3 coordinates of vertex 3
		* @param x4 coordinates of vertex 4
		* @param v1 value of vertex 1
		* @param v2 value of vertex 2
		* @param v3 value of vertex 3
		* @param v4 value of vertex 4
		*/
		Tetra(const float p1[3], const float p2[3], const float p3[3], const float p4[3], 
			float v1, float v2, float v3, float v4);

		/// Constructor
		Tetra(Vector3 coords[4], float vals[4]);

		/// Destructor
		~Tetra() {}

		static const int vert_neighboring[15][15]; 
		static const int six_simplex_neighbors[14][3];
		static const int six_tetra_index[6][4][3];
		static const int six_tetra_neighbors[6][4][4];
		static const int six_tetra_verts[6][4];

		/**
		* Compute the area of the isosurface (isoval=fx) within the
		* given tetrahedron
		* @param fx
		* 
		* @return float
		*/
		double isosurfAera(float fx);

		/// Extract isosurface within this tetrahedron
		void extractSurface(int id[4], float isoval, GeometryData* surface, EdgeCache& cache);

		///
		double getVolume() const { return volume; }

		/**
		* Compute the volume within a tetrahedron that has value <= fx
		* 
		* @param fx
		* 
		* @return float
		*/
		double innerVolume(float fx);

		/** 
		 * Compute the volume in a tetrahedron with value > fx
		 */
		double outerVolume(float fx)
		{
			return (volume - innerVolume(fx));	
		}

		/**
		* Compute the integral of multiple functions on the volume within a tetrahedron that has value <= fx
		* 
		* @param f1 function on vertex 1
		* @param f2 function on vertex 2
		* @param f3 function on vertex 3
		* @param f4 function on vertex 4
		* @param fx isovalue
		* 
		* @return float
		*/
		double funcIntegral(float f1, float f2, float f3, float f4, float fx);

		/// Compute the difference of functional interals between two isovalues
		double funcIntegralDiff(float f1, float f2, float f3, float f4, float isoHigh, float isoLow)
		{
			return (funcIntegral(f1, f2, f3, f4, isoHigh) - funcIntegral(f1, f2, f3, f4, isoLow));
		}

		/**
		* 
		* @param fv1 array of functions on vertex 1
		* @param fv2 array of functions on vertex 2
		* @param fv3 array of functions on vertex 3
		* @param fv4 array of functions on vertex 4
		* @param out output array
		* @param n number of functions
		* @param fx isovalue
		*/
		void multiFuncsIntegral(const float *fv1, const float *fv2, const float *fv3, 
			const float *fv4, float *out, int n, float fx);

	private:
		/// cooridnates of the 4 vertices
		/**
		* use double for better precision
		*/
		double x1[3], x2[3], x3[3], x4[3];
		/// values of the 4 vertices
		double v1, v2, v3, v4;
		/// original order of vertex ids
		int vids[4];

		/// variables used for tetra integral
		double area1, area2, midarea, volume, factor;

		/**
		* Sort the tetrahedron vertices such that v1 < v2 < v3 < v4
		*/
		void sortVerts();

		/**
		* Compute intermediate variables used by tetra integral
		*
		* called in the constructor
		*/
		void computeAreas();

		/** 
		* Compute the intersection point on a edge
		*
		* @return the index of the intersection point
		*/
		int interpEdge(int vid1, int vid2, double coord1[3], double coord2[3], double val1, double val2,
			float isoval, GeometryData* surface, EdgeCache& cache);

		Vector3 faceNormal(GeometryData *surface, int p1, int p2, int p3)
		{
			Vector3 v1 = surface->getVertex(p1);
			Vector3 v2 = surface->getVertex(p2);
			Vector3 v3 = surface->getVertex(p3);

			Vector3 t1 = v2 - v1;
			Vector3 t2 = v3 - v1;

			return t1.crossProduct(t2);
		}

	};

	// A pair of float and integer. Used for qsort
	class FloatIntPair
	{
	public:
		float   val;
		int     id;
		// Constructor 
		FloatIntPair(float _val = 0.0f, int _id = 0)
			: val(_val), id(_id){}

		// compare function used by qsort
		inline static int compare(const void* p1, const void* p2) {
			FloatIntPair* fip1 = (FloatIntPair*) p1;
			FloatIntPair* fip2 = (FloatIntPair*) p2;
			if (fip1->val < fip2->val) return -1;
			else if (fip1->val > fip2->val)	return 1;
			else if (fip1->id < fip2->id) return -1;
			return 0;	// same pair here
		}
	};

	// A pair of double and integer. Used for qsort
	class DoubleIntPair
	{
	public:
		double   val;
		int     id;
		// Constructor 
		DoubleIntPair(double _val = 0.0, int _id = 0)
			: val(_val), id(_id){}

		// compare function used by qsort
		inline static int compare(const void* p1, const void* p2) {
			DoubleIntPair* fip1 = (DoubleIntPair*) p1;
			DoubleIntPair* fip2 = (DoubleIntPair*) p2;
			if (fip1->val < fip2->val) return -1;
			else if (fip1->val > fip2->val)	return 1;
			else if (fip1->id < fip2->id) return -1;
			return 0;	// same pair here
		}
	};
}
#endif